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Abstract. Numerical evaluation of functional integrals usually involves a finite (L- 
slice) discretization of the imaginary-time axis. In the auxiliary-field method, the 
L-slice approximant to the density matrix can be evaluated as a function of inverse 
temperature at any finite L as Pt{0) — [f>i(P / L)] L , if the density matrix P\{j3) in 
the static approximation is known. We investigate the convergence of the partition 
function Zl(0) = Tr/3z,(/3), the internal energy and the density of states gh(E) (the 
inverse Laplace transform of Zl), as L — > oc. For the simple harmonic oscillator, 
gh{E) is a normalized truncated Fourier series for the exact density of states. When 
the auxiliary-field approach is applied to spin systems, approximants to the density 
of states and heat capacity can be negative. Approximants to the density matrix for 
a spin-1/2 dimer are found in closed form for all L by appending a self-interaction 
to the divergent Gaussian integral and analytically continuing to zero self-interaction. 
Because of this continuation, the coefficient of the singlet projector in the approximate 
density matrix can be negative. For a spin dimer, Zl is an even function of the coupling 
constant for L < 3: ferromagnetic and antiferromagnetic coupling can be distinguished 
only for L > 3, where a Berry phase appears in the functional integral. At any non-zero 
temperature, the exact partition function is recovered as L — > oo. 
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1. Introduction 

Functional integration is a long-established technique in quantum mechanics More 
recently, advances in computing power have allowed direct Monte Carlo evaluation 
of such integrals for many-body systems 0. Such algorithms are often based on an 
auxiliary-field functional integral, which is used in areas as diverse as strongly correlated 
electron systems |§, spin systems || and nuclear structure ||. The statistical mechanics 
of a many-body system on a d- dimensional lattice is mapped onto that of a classical field 
u(r, r) (the auxiliary field) in a d + 1-dimensional slab of extension < r < (3 = 1/kT 
in the imaginary time dimension. The many-body system reduces to a system of 
non-interacting particles moving in a time- dependent auxiliary field. To evaluate the 
integral over all time evolutions, it is necessary to sample the field at a finite number 
L of imaginary times — not necessarily uniformly or deterministically spaced — and 
extrapolate to the continuum limit L — > oo. The L — 1, or static, approximation maps 
the system onto classical statistical mechanics in d dimensions. The ground state in 
this approximation for many-fermion systems is a single Slater determinant, typically 
corresponding to the Hartree-Fock solution; for spin models it is the mean-field ground 
state. For L > 2, closed paths may enclose an area, breaking time-reversal invariance 
and thereby contributing a sign or Berry phase factor to the integral. This factor has 
some important consequences. It restores quantization: correlation between phases on 
neighbouring sites discriminates between the classically equivalent ferromagnets and 
unfrustrated antiferromagnets ||. The large- L limit must also restore symmetry if the 
auxiliary fields do not have the full local symmetry (such as in the Ising decomposition 
of the Hubbard model) 0. On the other hand, the resulting rapid oscillation of the 
integrand (the notorious sign problem) seriously restricts convergence of Monte Carlo 
simulations at low temperatures. The present author has shown how the distributions of 
the auxiliary fields tend to the appropriate quantum distribution (the Wigner function) 
with increasing L, while numerical convergence becomes increasingly problematical ||. 
For repulsive interactions, an imaginary auxiliary field is required, resulting in a sign 
problem even in the static approximation. 

Since numerical studies of the auxiliary-field functional integral are frequently 
hampered by the sign problem, it is of value to investigate toy models in which the 
finite-L approximants may be evaluated in closed form. The present work is a framework 
for discussion of these approximants, specifically for simple spin systems. This differs 
from finite-size scaling in real space; while a lattice truncated in real space is a cluster, 
and therefore physically realizable, the time-discretized system may possess unphysical 
properties vanishing only in the continuum limit. Indeed, in a number of examples 
the static approximants to the heat capacity and density of states are not positive- 
definite 0. The static approximation (and other finite-L approximations) give a saddle- 
point approximation, usually a variational overestimate, of the ground state energy, 
but are correct in the high-temperature limit. The heat capacity shows competition 
between the recovery of quantum fluctuations, which give a negative contribution at 
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low temperatures, and true thermal fluctuations, which give a positive contribution 
(exponentially small if there is a gap). 

To motivate this work, at this point we recall the path integral of a simple harmonic 
oscillator in the frequency domain (which does not suffer from the above problem). The 
partition function is [|TJ 
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We impose a frequency cutoff, restricting the function space to paths with L Matsubara 
frequencies (L odd): 
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The resulting Lth approximant to the partition function is [JTOj, |TT 
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with poles at (3 = 2irin/huj, (1 — L)/2 < n < (L — l)/2. The inverse Laplace transform 
of Zl gives the Lth approximant to the density of states: 

2^([(L-l)/2]!) 2 
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where is the Heaviside step function. The approximants have the following limits: 

9l (E) = B(E) (5) 

oo 

lim g L (E) -> V^-(n+ l/2)hu) (6) 

L ^°° n=0 

g L {E) ~E L -\E^0 + . (7) 

In particular, (|6|) verifies the emergence of the correct density of states in the continuum 
limit. Figure [I] shows the convergence of the internal energy Ul, obtained from Zl 
(0), to the exact result (hu/2)coth(hu /2kT) for any fixed positive temperature. The 
ground state energy vanishes for all L, and the low-temperature heat capacity is 
Lk; the zero-point energy is recovered with increasing temperature. In this case the 
finite-L approximants to the partition function represent physically realizable systems 
(ensembles of harmonic oscillators of the same frequency with a distribution of energy 
shifts). This is to be compared with the results to be shown in figure ^| and figure [], 
which do not exhibit this behaviour. 

Of more relevance to the present paper would be a time discretization of the path 
integral (|I]). The resulting approximants to the partition function have a similar form 
to (H), although the poles are non-uniformly spaced [1C, 12, H|. The approximant to 
the density of states has a less transparent form than (|4]), being quasiperiodic rather 
than periodic, but still converges with L in any finite energy interval. 




Figure 1. Internal energy for simple harmonic oscillator showing approximants with 
L = 1, 3, 5, 7, 9, 11, 13 Matsubara frequencies converging to the exact (bold) energy. 



Section ^| presents the auxiliary field formalism used in this work. Two case 
studies of toy spin models in section show how truncation of the functional integral 
gives a sequence of approximants which, although convergent onto the correct value, 
do not themselves represent any physical system. Section [| discusses possible wider 
applicability of the features of these models. 

2. Theory 

The general Hamiltonian with two-body interactions is of the form 

N N 

h = -J2 (^4 + K A D - E J^K A »- ( 8 ) 

Here A = {A^, \i = 1 . . . N} are single-particle operators generating a closed 
algebra; they will be spin operators in the examples studied here, but might, for 
example, represent hopping or pairing operators, cLc^ or cL-cL. To avoid notational 
complications, we assume the operators to be Hermitian. This can be achieved by 
changing the basis to A} + and i(A^ — A^). To obtain the functional integral, we 
separate the density matrix into L time slices, 

and apply the Hubbard-Stratonovich transformation to each time slice: 

eM-m/L) = • J ;'"^ Sg + g> ■ A ^ + Q(L-»). (10) 

/det(4jrJL//3) 



This is a formal expression, convergent only for a positive interaction matrix J. In 
general, one needs to reduce the matrix into positive, zero and negative blocks and treat 
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each separately, omitting the auxiliary fields in the zero block and using an imaginary 
auxiliary field in the negative block ||14j| . An alternative is to add a multiple of a 
positive matrix to J. In the latter case, this addition may correspond to a constant or 
one-body term, which can be absorbed into K, at the cost of introducing a fictitious 
self-interaction; the functional integration will have to work harder to remove this self- 
interaction. The coefficient of the additional term may be analytically continued or 



extrapolated to zero |15i ; it is this approach we shall use here. The Lth approximant 
to the density matrix is then 

Jd N u n exp(-/3u n ■ r l u n /AL) exp((3(K + u n ) ■ A/L) 



n=l 



(11) 



/ [det(4 7 rLJ//5)] 

This can be obtained from the density matrix in the static approximation 

PM = lPi(P/L)] L , (12) 

and tends to the exact density matrix as L —>■ oo. We shall subsequently refer to Lth 
approximant of "function" (obtained by replacing the exact density matrix with its 
approximant ( ]TT| ) ) as the L- "function" . 
The L-partition function is 

Z L (P) = Trp L (P). (13) 

Approximants to the internal energy may be computed directly from approximants to 
the partition function, 



which is an average of the one-body (auxiliary-field) Hamiltonian: 
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(14) 



(15) 



This will therefore tend to the expectation of the auxiliary field Hamiltonian at low 
temperatures, typically a mean field energy. 

Since it is difficult to extract the partition function from importance-sampled Monte 
Carlo calculations, the form flT4p is impractical. It is possible to calculate the energy as 
a thermal average of the true Hamiltonian 

U L ((3) = Z L ((3)- 1 Tr[p L ((3)H}. (16) 

This is a variational approximation to the ground state energy, which might be expected 



to be bounded below by the ground state energy. The example in section |37^ below shows 
that this natural assumption is not always justified for these approximants. The forms 
(H) and flTBj ) are not equivalent; the latter is usually a better approximation. Heat 
capacities will be defined as temperature derivatives of these energies, although these 
may be calculated in other ways |16| . 
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The partition function is the Laplace transform of the density of states gi{E). The 
L-density of states, gi,(E), is defined implicitly by 

/oo 
e-? E g L (E)dE. (17) 
-oo 

The spectrum is bounded below but, if necessary, the origin of E can be shifted to ensure 
that gi(E) = for E < 0. Such densities of states have been studied in the nuclear 
shell model JHj], although in that case the partition function is derived by integration of 
the measured energy ([16]) in (|TJ) and the inverse Laplace transform is computed within 
the saddle-point approximation (which is appropriate for a large density of states). The 
propagator may also be determined by a similar inverse transform of the density matrix. 
The L-partition function Z L (j3) converges pointwise to the partition function Z(fl) as 
L — > oo at any non-zero temperature. The L-density of states converges to the true 
density of states in the distributional sense: for any sufficiently smooth function /, 

/OO /"OO 
f(E)g L (E)dE = / f(E)g(E)dE. (18) 
-oo J —oo 

There are now two possibilities. If f(E)gi(E)dE is positive for all positive test 
functions f(E), the heat capacity is non-negative at all temperatures and we say 
that the approximant is physical; there can exist an Hermitian Hamiltonian with that 
thermodynamics. This is evidently the case for the harmonic oscillator discussed in 
section [I], although this is not related to the auxiliary-field functional integral. If the 
density of states is non-positive, then we say the approximant is unphysical. The spin 
models in the next section provide examples. 

3. Examples 

3.1. Single spin 

A single spin s with self-interaction, 

H = -JS-S, (19) 

may seem a trivial case, although a similar situation would arise in the study of a 
Hubbard model with degenerate bands and strong Hund's rule coupling. The static 
approximation to this has been discussed earlier ||. Although a scalar auxiliary field 
does not suffer from this problem, it violates rotational invariance |18[ . Clearly the 



exact partition function, internal energy and density of states are 

Z(P) = (2s + 1) exp(/5 Js(s + 1)) (20) 
U(/3) = -Js(s + 1) (21) 
g(E) = (2s + l)S{E + Js(s + 1)). (22) 
Applying the Hubbard- Stratonovich transformation gives the L = 1 density matrix as 
p x = (i3/4nJ) 3/2 J d 3 u exp(-(3u 2 /4J) exp((3u ■ S). (23) 
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Figure 2. Internal energy for single spin 1/2 showing approximants with L = 1 (top 
curve) to 10 (bottom curve) converging to the exact (bold) energy. 



This is rotationally invariant, and therefore a multiple of the unit matrix. The partition 
function in the static approximation follows from taking the trace of the exponential 
and performing the Gaussian integrals ||, giving 



Pi(/3) 



2s + 1 2s + 
The L-partition function is then 



- £ (l + 2m 2 /3J)e 



m 2 f3.J 



Z L {(3) = {2s + 1) (^-^ £ (1 + 2m 2 f3J/L)e m2 M L ^ 
this tends to the correct limit (|20"D as L —>■ oo. The L-energy ( 14]) is 

um - 



(24) 



(25) 



(26) 



E^=- S (l + 2m 2 /3J/L)e m2 ^/ L 
This is a monotonically decreasing function of temperature, falling from the saddle-point 
value of — Js 2 at T = to the correct value of — Js(s + 1) at high temperatures. Figure ^| 
shows the energy for spin 1/2. The thermal average of the Hamiltonian ( |T6"D is trivially 
Ul{(3) = —Js(s + 1) at all temperatures. 

This negative heat capacity implies a non-physical density of states. For spin 1/2, 
the approximants to the partition function are 

Z L {(3) = 2{l + (3J/2L) L ^ J l\ (27) 

converging to 2e 3/3J/4 as L -> oo. Its inverse Laplace transform is the L-density of states 

L 

(28) 



J d 



which involves L derivatives of the delta function at the saddle-point energy. The formal 
limit 



lim g L (E) = 2exp ( 1 + ^ ) 6{E + J/4) = 25(E + 3 J/4) 
L— >oo \ z ah, / 



(29) 
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is correct when applied to a sufficiently good test function (such as a polynomial). In this 
way expectation values are correct to 0(L _1 ), even though the only spectral point is at 
the classical rather than the quantum ground state energy, and is in error by 0(1). For 
larger spins, this single singularity becomes a discrete spectrum of singularities between 
E = —Js 2 and E = — J/4 (half-odd- integer spin) or E = (integer spin), all higher 
than the true energy. 

There is a suggestive but probably fortuitous resemblance between the partition 
function for spin 1/2 and the g-Laplace transform, recently-introduced in the context of 
non-extensive statistical mechanics . One version of this g-Laplace transform defines 
a g-partition function as 

Z q (/3) = / g(E)[l + (1 - ^fif/NdE, (30) 

where q corresponds to 1 — 1/L. Negative heat capacities are found in this theory [p0[] 
(although the correspondence between Z q and thermodynamic potentials differs from 
that in standard thermodynamics). 



3.2. Spin 1/2 dimer 



The highly non-physical behaviour of the approximants to the heat capacity above stems 
from the emergence of the quantum fluctuations (a negative energy contribution) with 
increasing temperature. For a single spin in zero field there are no compensating thermal 
fluctuations. We therefore investigate the spin- 1/2 dimer, 

H = -J'{S X -S x + S 2 - S 2 ) - 2JS X ■ S 2 . (31) 

The self-interaction is added to ensure convergence. The integral ([10[) only converges 
for J' > |J|, but is analytic in the matrix elements, allowing continuation to J' = 0. 
Manipulation of the Gaussian integrals eventually gives the L-density matrix as 

'5 _ j/V P(J+JT \ BJ/2L ,(1 ,J^_ _ P(J ~ J"' 

6 6 J 2 3LJ I 1 6 6 J 2 



6J 2 
J 



a 



3LJ 

a (3J/2L 



-0J/2L 



6LJ 



+ 1 



J 



12 



1 L 



(32) 



2J 2 / \2 2J 2 2LJ 

This is an entire function of both J and J' . Taking J — > gives the direct product of 
two one-particle density matrices (|25[) . More importantly, we can remove the interaction 
by setting J' = to obtain 

'\_PJ_ 
6~6L 



5 pr 

6 3L 



e /3J/2L + 



-0J/2L 



Pi 



+ 



± e PJ/2L 



3 _ PJ \ -pj/2L 

2 2L 



(33) 



where P\ and Po are projections onto the triplet and singlet subspace respectively. In 
the large-L limit we recover the correct density matrix: 



lim p L (p) 

L—>oo 



lim 

L— >oo 



1 + 



2L 



Pi 



1 



3PJ 
~2L 



Po 



(34) 
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= e PJ/2p i+e -3PJ/2p o (35) 

It is not possible to distinguish ferromagnetic from antiferromagnetic coupling in 
the thermodynamics for L < 3, where the paths do not enclose an area. Thus the L- 
partition function is an even function of J for L — 1, 2. A high-temperature expansion 
(Maple) verifies this and shows that the second moment of the density of states is correct 
for all L: 

ZM - 4+— ^ + ^ +■ • • -( 36 ) 

Figure ^ shows approximants to the internal energy for ferromagnetic and 
antiferromagnetic coupling. Ul, as calculated from the partition function (|T4|), is always 
equal to its mean field value —\J\/2 at T = 0. This is the correct energy only for the 
ferromagnet. U\ and Ui, as already discussed, cannot distinguish the ferromagnet and 
antiferromagnet. The L-heat capacity is negative at low temperatures and positive at 
higher temperatures. Ul, as calculated from the thermal average of the Hamiltonian (Pj) 
is, as expected, a better approximation than Ul for L > 1, giving improved estimates 
of the antiferromagnetic ground state energy, although it still shows a small region of 
negative L-heat capacity. 

One at first surprising feature is that in the ferromagnet Ul falls below its variational 
bound —J/2 for odd L. This is due to the unphysical form of the L-density matrix itself, 
and not just to its temperature dependence. As a result of the analytic continuation 
to J' = 0, the coefficient of the singlet projector in (|33"D is negative at low temperature 
for odd L, representing a negative weight for the singlet state. The coefficients of both 
projectors are always positive when J' > \J\, the parameter region for which the integral 
(|10|) converges. Direct application of the Hubbard-Stratonovich transformation to the 
Hamiltonian for J' = would require an imaginary field coupled to Si — S2, leading to 
a similar non-classical weight. 

The L-density of states is again non-physical and is symmetric for L < 3 for the 
reasons discussed above; there are L + 1 singularities in —J/2 < E < J/2, involving L 
derivatives of the delta function, for example 

g x {E) = 2 (5(E + J/2) + 5(E - J/2))+J (5'(E + J/2) - 5'(E - J/2)) .(37) 

For large L, we obtain the correct result (the limit to be understood in the distributional 

sense) 

lim g L (E) = 35(E + J/2) + 5(E - 3J/2). (38) 

L^oo 

4. Discussion 

In all the above, the approximants to thermodynamic functions have error 0(L _1 ) 
at any fixed non-zero temperature. In practice, more careful Trotter decompositions 
and truncations of the density matrix may accelerate the convergence in Monte Carlo 



simulations JL(| [21], p2| . The models discussed, being analytically soluble for all 



discretizations, are not representative of real applications but may provide a useful 
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Figure 3. Lth approximant to the energy for spin-1/2 dimer with (a) ferromagnetic 
and (b) antiferromagnetic coupling. The bold line is the exact energy, the full lines are 
Ul, the derivative of InZ^ (|l4|), and the dashed lines are Ul, the thermal average of 
the Hamiltonian (|l6|). Curves shown for L = 1 . . .4, annotated by L. The energies U\ 
and Ui are coincident. 

test of methods. The main outcome of this work is a pointer to possible difficulties in 
the use of finite discretizations: underestimated (or negative) heat capacities and non- 
physical spectral functions. The effects might be largest in strongly correlated systems, 
or systems with an excitation gap, where the lowest auxiliary field state is a poor 
approximation to the true ground state. In that case the density of states must suffer 
substantial distortion to provide the correct thermodynamics; the true ground state 
energy lies outside the approximate spectrum. More accurate energies are obtained 
from the thermal average of the Hamiltonian than from the derivative of the partition 
function. 

This work has been in some sense complementary to the sign problem, which 
can arise for L > 3 or for repulsive interactions fL4}| : the weight in the functional 
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integral (|1(J) need not be positive, although physical quantities are correctly obtained. 
This will arise if the correlations to be calculated are incompatible with a positive 
distribution for the auxiliary fields ||. In this case low-L approximations lead to 
unphysical results, characterized by non-positive distributions in the energy domain. 
The examples discussed are those with the most acute sign problem: the auxiliary field 
couples to operators (such as spin components) whose equal-time commutators do not 
vanish. 
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